{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import string\n",
    "#https://stackoverflow.com/questions/19726663/how-to-save-the-pandas-dataframe-series-data-as-a-figure\n",
    "import six"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "tbpath = \"../../fits/\"\n",
    "productpath = \"../../postfit_derivatives/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "models = [\"fulllinearmodel_fit_table.csv\",\"reducedlinearmodelNegBinom_fit_table.csv\",\n",
    "          \"reducedlinearmodelq0_fit_table.csv\",\"reducedlinearmodelq0ctime_fit_table.csv\",\n",
    "         \"nonlinearmodelq0ctime_fit_table.csv\",\"nonlinearmodel_fit_table.csv\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "rois = []\n",
    "for model in models:\n",
    "    df = pd.read_csv(tbpath + model) #get rois in all tables (some may have failed)\n",
    "    rois += list(df.roi.unique())\n",
    "  \n",
    "rois = list(set(rois))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gambia\n"
     ]
    }
   ],
   "source": [
    "dfreport = pd.DataFrame(columns=['Region', 'Model','median CAR0','CI lower','CI upper','median CARc','CI lower','CI upper','delta weeks'])\n",
    "k = -1\n",
    "for roi in rois:\n",
    "    k += 1\n",
    "    model = \"fulllinearmodel_fit_table.csv\" #dfwaic_sorted.loc[dfwaic_sorted.roi==roi,'model'].values[0]\n",
    "    df = pd.read_csv(tbpath + model)\n",
    "    theta = \"car (week 0)\"\n",
    "    measure1 = df.loc[(df.roi==roi)&(df['quantile']==0.5),theta].values[0]\n",
    "    measure2 = df.loc[(df.roi==roi)&(df['quantile']==0.025),theta].values[0]\n",
    "    measure3 = df.loc[(df.roi==roi)&(df['quantile']==0.975),theta].values[0]\n",
    "    x = [roi,model.split('_fit_table.csv')[0],np.round(measure1,4),np.round(measure2,4),np.round(measure3,4)]\n",
    "    #find latest week with data\n",
    "#     maxweek = 11\n",
    "    for i in np.arange(11,0,-1):\n",
    "        theta = 'car (week '+str(i)+')'\n",
    "        x2 = df.loc[(df.roi==roi)&(df['quantile']==0.5),theta].values[0]\n",
    "        if np.isfinite(x2):\n",
    "#             print(x2)\n",
    "            measure1 = df.loc[(df.roi==roi)&(df['quantile']==0.5),theta].values[0]\n",
    "            measure2 = df.loc[(df.roi==roi)&(df['quantile']==0.025),theta].values[0]\n",
    "            measure3 = df.loc[(df.roi==roi)&(df['quantile']==0.975),theta].values[0]\n",
    "            x += [np.round(measure1,4),np.round(measure2,4),np.round(measure3,4),i]\n",
    "            break\n",
    "    try:\n",
    "        dfreport.loc[k] = x\n",
    "    except:\n",
    "        print(roi)\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0864\n",
      "0.025    0.03568\n",
      "0.500    0.08640\n",
      "0.975    0.19412\n",
      "Name: median CAR0, dtype: float64\n",
      "0.186\n",
      "0.025    0.01300\n",
      "0.500    0.18600\n",
      "0.975    0.43848\n",
      "Name: median CARc, dtype: float64\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2oAAAFCCAYAAACASy55AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeVxV1f7/8ddCFARnwXnA0JyuSTmbIuKAzWpqg2XkUHrVLKeyTLGyvFpXcyibFMvqp1+nbqZmili3nLW6Gdqg6E3NKWdxwv37Yx/O5cABjniQg76fjwcPYO211/7sxeGc8zlr7bWNZVmIiIiIiIiI7/DL7wBERERERETElRI1ERERERERH6NETURERERExMcoURMREREREfExStRERERERER8jBI1ERERERERH6NETUSkgDHG1DDGLDHGHDbGWMaY+Fy04fF+xpg4R/2wKz1OujZiHW1EZVfmKHd7fsaYIGPMVGPMXmNMqjEmObfxiIiI+Dr//A5ARESuWDxwCzAe+BP4PV+j8b543J/fs8Bg4HXgR+BUfgQnIiJyLShRExEpQIwxAUBrYLplWa/ndzxX6SPg/wEX0gpyOL8OwH8syxpx7UIUERHJH5r6KCJSsJQHDPBXfgdytSzLSrUs65xlWZfTFWd3fhWyKL8qxpji3m5TRETkailRExEpIBzXau1x/DrWcf2W8xovY4y/MeZZY8zPxphzxpijxpjFxpgGHrbvZ4wZZYzZ7dj/J2NMz1zE2c8Ys8MYc94Y85sx5mns5CtjvdgM8Wd1frHGGAuoAbRJVx6Xrq3GjnM94jjuTmPMC8YY/wzHTDTGJBtjbjLGLDDG/AWcTLfdGGMGGGO2GGPOGmNOG2PWGGPaZmgnLC0GY8zdxphNjj47YIyZlPG4jn1qGmNmG2P+MMZcMMbsN8Z8ZoxplKGeR+ciIiLXNz3pi4gUHO8A3wOTgcXAIkd5kuP7x0AP4CvgbewRqIHAOmNMa8uytuXQ/j+BIcDXjmOUA2YAuzwN0JGUTQZ+AJ4HgoDhwCEPds/q/H4FHnWUH8G+dg3s69QwxtzlqPsb8Ab2qFsL4CUgAuie4TjFgLXAt8ALjvNM8xHwELAAmA0EAD2Br4wxXS3L+leGtu4E/g7MBGYB9znO9xjwarp+aQysBgoDHwA/AWWANkBLYEsuz0VERK5TxrKs/I5BREQ85Fh5cTcwzrKsuHTlHYCVwHzgQcvx5G6MaYidBKyzLKt1uvoWMMeyrFjH77WxE741QEfLslId5bcBm7FHxGpYlpWcTWylgH3Yo2KNLcs66yivAuwAgoG2lmUlOspjsZOh9GVuz8+xLRlItiwrKl1ZIJAM/AJEW5Z1Kd22Z7CTz/TtJ2InR+Mtyxqdof0u2EnSk5ZlvZuu3B9YD5QFbrIsy0oX51mgflq/GGMM8B+grGVZFTOU1QSaWpb1Y4bj+lmWdflKz0VERK5vmvooInJ96OL4Pt5K9wmcZVk/AJ8DrYwxodnsfx92MvbPtCTNsf9W7BE6T3TEHkGbkZakOdr4A3u0Ly90wL6ubTZQyhgTkvYFLEsXV0buFmJ5BHslySUZ2imF3YdhQK0M+yxJn7w6+n4NUMEYU8xRHAHUB2ZnTNIc+6Rdo5fbcxERkeuQpj6KiFwfagCX+d80yPS2A50ddQ5nsf9Nju873Gz7Gc8ShJzayAt1Hd9nZVOnfIbfD1uWdTyLtooDB3No65d0v7ubFnrU8b0scJr/JXc5TT3NzbmIiMh1SomaiIgUZGmLlIzAvr7Nnf0Zfj/rtpbd1mHg4WyO91OG31Pd1nKNzVO5ORcREblOKVETEbk+7MKezl4XxyIb6dRzfN+dw/4Adch8A+16eCZ9G6tz2caV+tXx/YxlWau80NbNwHrLsk5fZVvppY3ARXhwfPDOuYiISAGna9RERK4PSxzfRzkWrwDAGPM34F7g35ZlZTXtEeBfgAUMNcYUSrf/bUB7D2P4CkgBBhpjgtK1UYXsR6muxpfYK0o+Z4wpk3GjMaao8fw+aR9ivy6+5m6jMSa30w5/wJ5+2tsYU99Nu2l/L2+ei4iIFHAaURMRuQ5YlvWVMWY+8CBQ2hizlP8tz38OeCqH/XcYY2YAg4AEY8xC7GXrB2EnGrd6EMMxY8yL2At1fGeM+RB7cZH+2KNFObZxpSzLOmOM6YWdqO40xszCXtq+FPbIXlfshVYSPWhrgTFmNjDIkaAuxb4dQBXsJfJr8r/r8K4kRssY8zj2KONGY0za8vylsFegXAFM8+a5iIhIwadETUTk+tET2ArEYt+D6wz2/cJetCzrPx7sPwT4E3gCmISdXA3EXgzDoyTLsqw3jDGngaHYI1P/xU7cTpD9Ihm5ZlnWl8aYJsBz2Cs3hmLfx+x37CXtM620mE1bvY0xa7D7YBRQBLtPtjp+z22Mmxwxvoh9r7v+2EngRuz7uXn9XEREpGDTfdRERERERER8jK5RExERERER8TFK1ERERERERHyMEjUREREREREfo0RNRERERETEx/hMomaMudkY85IxZr0x5rAx5pQx5ntjzAvGmGA39WsbY5YYY44ZY84YY74xxkRf4TFz1YYxZpQx5v+MMbuMMZYxJjmH+s2MMasc53TSGLPCGJPTjU+93kZB4uhXd18e34TWGHOnMeY7x9/2L8ffrEZexu2LHI/zj40xScaYE8aYs8aYHcaYfxpjKnqwf2FjzExjzBZjzBFjzHljzG5jzDxjjNeXW/dlV/o8lUUbidk8vhvn9Tn4Em/0p6Md/a87XOnrk5v9Y7N5fKZ9Vc6j8H3K1falm/b8jDHrHG0t9VKYBYYXHpuVHW2sNcYccPy/bzfGTDLGlM2jsH2WNx6fxpgnHe8PdhhjUo0xN+wKg976fzfG9DLGbDPGpBhjDhpj3jfGhOY6Ll9Z9dEYMwF7Geh/AeuBi0Bb7GWMfwSaW5aV4qgbjr2k8SVgCvayz/2AvwF3WJa1yoPj5boNxwP5L+zlmhsBJy3LCsuibnPse97sA6Y7igdh35+opSdLZnujjYLG0cffAO9m2HTRsqx5HuzfFViAff+n94CSwNNAKtDYsqz93o3Ydxlj2gEvYP9f/YH9mG8APA6cBCIsyzqUzf7B2Eu8fwfsAk4B1Rz7VwA6WZaVkJfn4Cuu5HkqmzYSgfrAM242L7Ms6y9vxuzLvNSf+l9P50pen7LY/yagpZtNFYGJwDbLsm7zQqg+72r70k17g7BvWVEM+MKyrLu9EWdB4YXHZn/gTeAL4N/Yr0VNsW9H8ifQxLKsP70bte/yxuPTkYyUBbYBNYAqlmUZ70ZaMHipP5/Bvo3KWuAT7HtwDgX2AE0tyzpzxYFZluUTX0BjoKSb8lcACxiUrmw+9otwRLqyYo6O2IkjAc3heLluA7gp3c8/AcnZ1N2I/Wa4crqyyo6ylR72zVW3UdC+HH/z+FzuWxg7qd0DFEtXHuH4m7+b3+fnC19Ad0c/j8zl/hWx31gvy+9zuYZ95vHzVDZtJGb3nHEjfV1tf+p/3W2fePz6dIXtjnL8TQbm9zkWxL7EfsN2EvtNmwUsze/zK2j9if0BVwU35X0dffp6fp9jQepPx35hgJ/j56WAld/nVVD7EwjBvnfpRqBQuvJ7HI/P53MTl89MfbQsa7NlWSfcbEobPfkbOD/dvxdItCzr+3T7nwbeB24GmmR3rKttw7KsXZ6ckzGmpqOd/7Msa1+6/fcB/we0N8ZUyOs2CjJjTBFjTLEr3K0NUAl43/E3BcDxt04EHjDGFPZelAXWHsf30rnc/xBw7ir2L3A8fZ7yhGMaVAljzA356SV4pT/1v56Bp69PV8LxGO0NpAAfe7t9X+XlvpyBPSPhTS+2WaBcbX9alrXdcj9idsXPv9cDbzw+LctKtizrsjfiKei80J+dgSBgmmVZqena/Rz7f/+R3DTqM4laNqo4vh90fL8FCADWuam73vE920TNS214Iq2NrI5jsIdX87qNgqobcBY4ZYw5ZIyZZowp6cF+OfVZCexk/IZijAk0xoQYY6oYYzoC7zg2LfNw/0KO/SsYY5pgD+sX83T/61zG56mcVAZOY0+5Pm2MWWSMqZMnkRVMnvan/tevjTZATWChZVnH8zuYgsYY0w37U/X+6d/Aiddc6fOvSF7I6fWoTi4GHvC/qpDymDGmEPAi9jU1nziKKzm+73OzS1pZThc6e6MNTxSkWH3NRuwRw9+w32zdiX1dXhtjTMv0n5674WmfbfdSrAVFX2Baut+TgUcsy/rGw/3rAumvhzyBfb3Fa16JroDK4nkqO7uBb7GvwUoFmmE/ttsZY1pZ1+E1p1fiCvtT/+vXRh/H9/fzNYoCyPHh4lTgHcuy1udUX3JlnOP7nHyNQm50Ob0eGUedX66kUZ9O1LAX+WiBPa9zp6MsyPH9vJv65zLUyYo32vBEQYrVp1iW1SxD0YfGmB+B8cAQx/es3JB95oElwA7sUbBbsaf/hlzB/ruBDkAR7E/XH8FeuCEA+031jcrd81SWLMt6PEPRAmPMv7Cn6v0Tu49vZFfSn/pfz2PGmFLA/cBvlmWtze94CqCJ2LOXRuV3INcjY8ww7Out37VukEWtxGflyeuRzyZqxpiXsT9lfteyrPSf2J91fA9ws1tghjpZ8UYbnihIsRYEk4CxwF1kn6ipz9ywLOsP7FUfAZYYYxYCm4wxQRn+x7La/wzgXA3VGDMLe3WkhUCnPAjZ52XzPHVFLMv6xhjzNdDWGFPUymGlw+tVLvpT/+t572GgKPBBfgdS0BhjWmOvJv2opox6nzGmL/b7gi+wnzdE8lP616OMr+G5fj3yyWvUjDFxwGhgNtA/w+a0pZbdTfdLK3M37OjtNjxRkGL1eZZlXcTuj5xGgdRnHrAs60fsJXn/nsv9TwOLgBjH7S5uKDk8T+VGMlCIG2hxlvRy2Z/6X897fbBHzOPzOY6CaDr2bSM2GGNqpn05tgU5fr+SWQ3iYIzpjX37npXA/Y73ByL5KafXIytdHY/5XKLmeLEeiz3XuK/lWNsynf9gDyu2cLN7c8f3zTkcxhtteGKT43tWx7GALdegjeuCMSYQ+6LhnC4YzqnPTnKFc4SvY0WBMle5P1fZRoHjwfNUbtTCfkN8w9xHLc1V9Kf+1/OQMSYCuA37nl83zP2pvKg69q0ifs3wBfb9An8F4vIlsgLMkaS9jz3Do7NlWe6mmolcazm9Hu3MYX0Ft3wqUTPGjMF+sf4I6O1uyVDHSX4ORBljGqbbtxj2Ygm/Yi9EkVZe0hhTJ/2nVlfaRm5ZlvUbdsLX3RiTdpEhjp+7AwnpX/wcK+rVSb+y4ZW2cT0wxpTNYtPL2NN1P09Xt6Kjz9LP+10LHAD6pl9hx/G3jsK+1cEN8+lbVrdvMMa0xV7OeH26skz9aYwJNcZkeq5wtNsde/XCG2axBk+epxz13PVlScdiGRnr3gXcDnxlWda5jNuvZ1fTn+h//aoYY6o5+jSrWxj0dXzXtMccZNGXvbCfIzN+gf0Ba3e0QItbWT02jTGx2De2TwDuu9GeL3PLg/91uQJZ9Odn2FMeB6V/nTfG3APcRC5vbWK880Hw1TPGDMSeJrAXe8WvjC/WBy3L+spRtyZ2InURmIz9qWk/oAFwl2VZX6ZrNxZ7Ks04y7Li0pV73IabWB/F/qQMYDD24gpvOH7fY1nWR+nqtgTWYF8bNC3dPuWB2y3L+iFd3TjsNyyPW5YVn5s2rgfGmMnYnz6swX48FMNe9bEtsAFom3YNjzEmHnjMUZaYro3u2PdW+QH7Sb0E8Az2CGSj9Peku94ZYxZj35w6AfveaYHYt3R4EHu+dFTa/QTd9acx5mngaWAx9oIiF7CXPH8Me5peX8uyZl27M8o/V/g8FU/mvuyMvWBI2n1VLgFNsRdm+Qv7//mGGQG62v50lOt/PZ0rfH1KxF56v4ZlWckZ2gnEToLPAtVuxGXlvdWXbtq1sEcp7/Z2zL7savvTGHMv9uvQSWAkma8DOm1Z1pK8it/XeOPx6Ugi0gYsHgFqYz8XAxy3LGt6XsXva7zUn8OA17EXB/sUe8rjMOC/QJPcjKjl+53A076w579b2XwlZqhfFzt7PY79QvJvoL2bdmMd+8e52eZRG272S/Q0Tkf9FsBq7JGHU8CXwG1u6sU52ojNbRvXwxdwn+P89mGvlHMG+B54HgjM4nET5aadu7FHi84Cx4AFQHh+n18+9GcPYCn2E8U57Be3HdhJf7Wc+hM7qZuLPdJ8GjtR+y/2m+OW+X1+17gvPX6eyqIv6wLzgd8dfXne8fMMoHJ+n19B68902/S//r++8Pj1KV3dMDftPOzYNj6/z6mg96Wbdi1gaX6fX0HrT/73Himrr+T8PseC1J+O8uyeg9WfuXvujMX+4PAccAiYBZTLbVw+M6ImIiIiIiIiNp+6Rk1ERERERESUqImIiIiIiPgcJWoiIiIiIiI+RomaiIiIiIiIj1GiJiIiIiIi4mOUqImIiIiIiPgYJWoiIiIiIiI+5rpM1IwxT+R3DDkpCDEWJOpP71J/epf603vUl96l/vQu9af3qC+9S/3pXdeqP6/LRA0oCA/GghBjQaL+9C71p3epP71Hfeld6k/vUn96j/rSu9Sf3qVETURERERE5EZkLMvyuHJISIgVFhaWd9F4yeHDhwkNDb1mx9t5dCcAtcvW9nifax1jQXL07FEAygaV9Xgf9ad3qT+9S/3pPepL71J/epf603vUl96l/vSuLVu2nLYsq3heH8f/SiqHhYWxefPmvIqlwIqKjwIgMTYxX+O4XsR/Hw9AbERsvsYhIiIiIpKRMWbntTiOpj6KiIiIiIj4GCVqIiIiIiIiPkaJmoiIiIiIiI9RoiYiIiIiIuJjlKiJiIiIiIj4GCVqIiIiIiIiPkaJmoiIiIiIiI+5ovuoiYiIyPXr3LlzHD58mHPnznHp0qX8DkdE5Jrw9/cnMDCQ0NBQAgMD8zscp+s2UXt/3YKr2r9vi25eikRERMT3nThxgoMHDxIaGkqFChXw9/fHGJPfYYmI5CnLsrh06RKnT59m7969lC9fnpIlS+Z3WMB1nKiJiIiI544cOUKVKlUICgrK71BERK4ZYwyFCxemdOnSBAQE8Oeff/pMoqZr1ERERIQLFy5QtGjR/A5DRCTfFC1alPPnz+d3GE5K1ERERARAUx1F5Ibma8+BStRERERERER8jBI1ERERERERH6NETURERCSfJCYmYowhPj4+2zKxGWOIjY3N7zBErgklaiIiIiJyw9i5cyedO3emdOnSBAcH07p1axISEvI7LJFMtDy/iIiIiA+JjIwkJSWFwoUL53co153ff/+dli1b4u/vz8iRIylZsiTvvfceMTExLF++nPbt2+d3iCJOStREREREfIifnx+BgYH5HcZ1adSoURw/fpwtW7YQEREBQK9evahfvz4DBw5kx44dPrfyn9y4NPVRREREbijx8fEYY1i9ejUvvfQS1atXp2jRojRr1oz169cDsHbtWlq1akVwcDAVK1bk5ZdfdtvW5s2b6dKlCyEhIQQEBFC7dm3Gjx/PpUuXMtX97LPPuPXWWwkMDKRq1aq8+OKLXLx4MVM9d9eoXb58mfHjxxMZGUmFChUoUqQI1apVY8CAARw9etRl/+TkZIwxxMXFsXTpUpo0aUJgYCAVK1ZkxIgRbmPLqG3btoSFhbmUffrppxhjaNiwoUv522+/jTGGDRs2OMssy+Ltt9+mUaNGBAUFUaxYMdq2bcuaNWvcHm/evHm0atWK4sWLExQURLNmzViwYEGOcQJs3bqVChUqUK9ePfbu3ZtlvTNnzvCvf/2LqKgoZ5IGUKxYMfr27csvv/zCpk2bPDqmyLWgRE1ERERuSM899xxLlixhyJAhjB07ll27dtGxY0eWLFlC165dad26Na+//jp16tRhzJgxzJ0712X/L774gttvv51ffvmFYcOGMXXqVFq0aMGYMWN46KGHXOouXryYLl26cOLECcaMGcPgwYNZsGABzz77rEexXrhwgUmTJlGrVi1GjBjB1KlT6dChAx988AFRUVFcuHAh0z7Lli2jd+/e3HHHHUyePJmGDRvy+uuvM3HixByPFx0dzZ49e/j999+dZatXr8bPz4///Oc/HDlyxFmekJBAiRIlaNy4sbPs0UcfZdCgQdSsWZOJEycybtw4Tpw4QYcOHfjXv/7lcqzRo0fz4IMPUrx4cV5++WUmTJhAUFAQ3bt3Z8aMGdnG+eWXX9KmTRvCw8P597//TbVq1bKs++OPP3L+/HlatGiRaVvz5s0BlKiJT9HURxEREbkhpaamsn79eooUKQJAvXr1uO++++jevTvr1q1zJh59+vShevXqzJgxg0ceeQSAc+fO0adPH5o1a0ZCQgL+/vZbqieffJKGDRsydOhQEhMTiYqKIjU1lSFDhlCmTBk2btxISEiIs+4tt9ziUawBAQEcOHCAokWLOsv69+9Py5Yt6du3L0uWLKFHjx4u+2zfvp3t27c7R8b69+9PgwYNmDZtGs8//3y2x4uOjmbMmDEkJCQQHh4O2AnZww8/zNy5c0lISKBHjx5YlkViYiKRkZEUKlQIsJPSjz/+mHfeeYcnnnjC2eaQIUNo3rw5Q4YM4Z577sEYw9atWxk/fjyjRo3i1VdfddZ96qmn6Ny5M6NGjaJXr14UL148U4wfffQRffr04c477+TTTz916Rt39u/fD0DlypUzbUsr27dvX7ZtiFxLStREREQkS0+veJrv//w+v8NwEVEhgimdplx1OwMGDHAmaQCtW7cGoFmzZi6jQ0WKFKFp06Z8++23zrKvvvqKgwcP8tprr3H8+HGXdu+8806GDh3KypUriYqKYsuWLfz3v/9l+PDhziQNoGTJkvTv3z/HpAnsZenTEpHU1FROnTrFpUuXiI6OBmDDhg2ZErXOnTu7TF80xtC2bVumT5/O6dOnKVasWJbHa9q0KcWKFSMhIYF+/fqxZ88edu/ezfTp0/n+++9ZvXo1PXr0cI6upcUBMHfuXIoXL07nzp1dRt4A7rnnHuLi4vj111+5+eab+fjjjzHG8Nhjj2Wqe++99/LZZ5+xbt06Onbs6LJtwoQJPP/88zzxxBPMmDHDmSRm5+zZs4Cd9GaUdk1gWh0RX6BETURERG5IN910k8vvpUuXBqBGjRqZ6pYuXdrlWrCkpCQAevfunWX7Bw8eBGDXrl0A1KlTJ1OdevXqeRzv/PnzeeONN9i2bVuma9uOHTuWqX7G8wMoW7YsAEePHs02UStcuDCtWrVyXlO2evVq/P39iYyMJDo6mmXLlgE4l7VPn6glJSVx6tQpypcvn2X7Bw8e5OabbyYpKQnLstz2Tfq66S1atIhTp07Rr18/Zs6cmeV+GQUFBQFw/vz5TNvOnTvnUkfEFyhRExERkSx5Y+TKV2U1CuPJ6IxlWQBMmjTJZWGK9CpVqpT74DJYtGgRDzzwAE2bNuXNN9+katWqBAYGkpqaSqdOnbh8+XKmfbI7j7T4sxMdHc2KFSvYvn07CQkJNGnShGLFihEdHc3UqVPZu3cvCQkJhISEuEzhtCyL0NBQPvnkkyzb/tvf/uasa4xh+fLlWcZbv359l9+bNm1KcnIyCxYs4IknnnAZ/cxO2t/D3fTGtDJ30yJF8osSNREREZErVKtWLQCCg4NzvPdW2sjWjh07Mm37+eefPTreRx99RGBgIGvWrHEZ9XHXprekjZKtXr2ahIQE+vTpA0BUVBSFChVi5cqVfP3113To0MFlSftatWrxyy+/0Lx582xH7dLqrlixgmrVqlG3bl2P4qpSpQpz5swhOjqa9u3bs2LFCudiINlp0KABAQEBrFu3LtO2tNU+PU36RK4FrfooIiIicoViYmIoV64cEyZM4K+//sq0PSUlhVOnTgHQqFEjqlSpwuzZs12uwzp58qTHU/cKFSqEMcZl5MyyLF555ZWrPJOs3XrrrZQuXZqZM2dy4MABZ+JWsmRJbrvtNiZPnsyJEydcpj2CfV+yy5cvM2rUKLftpp/K+OijjwLw/PPPk5qamm3d9CpXrszatWupVKkSHTt2dLl+MCvFihXjnnvuITExkR9++MFZfvr0ad5//31q1apF06ZNc2xH5FrRiJqIiIjIFQoODubDDz+kc+fO1K5dm969e1OzZk2OHz/Ojh07WLRoEYsXL3aOPk2ePJkePXrQtGlT+vXrh7+/P7NmzaJs2bLZ3vsrTbdu3Vi4cCHR0dH06tWLixcvsmTJkjxd/MLPz482bdqwZMkSAgMDadmypXNbdHQ0//jHP5w/Z4z18ccfZ/r06WzdupW7776bkJAQ/vjjD9atW8dvv/3mvG6vSZMmxMXFERcXR0REBN27d6dSpUocOHCALVu2sGzZMre3HgCoUKECiYmJtG/fnk6dOrF06VLatGmT7Tm99tprrF69mo4dO/LMM89QokQJ3nvvPfbt28cXX3yhm12LT9GImoiIiEguxMTEsGnTJmJiYpg7dy4DBw7k9ddfJykpiaFDh7pct9WtWzcWLFhAiRIliIuLY+rUqXTr1s2Z7OTkwQcf5N133+X06dMMHz6ciRMnUrt2bb788su8Oj3gf0lYy5YtXVZLbNeuHWCPbNWuXTvTfrNmzeLDDz/Ez8+P1157jcGDBzNnzhyKFSvGa6+95lJ37NixLF26lEqVKjFlyhQGDhzIu+++y/nz55k6dWq28ZUrV441a9ZQs2ZN7rzzTlavXp1t/Zo1a/Ltt9/SvHlzJkyYwPDhwwkODmbFihXExMR41Cci14rx5GLSNI0bN7Y2b96ch+F4z/vrPLubfVb6tujmcd2o+CgAEmMTr+qYYov/Ph6A2IjYfI1DRORGkpSU5PE1QiIi1ytPnguNMVssy8rzCxo1oiYiIiIiIuJjlKiJiIiIiIj4GCVqIiIiIiIiPkaJmoiIiIiIiI9RoiYiIiIiIuJjlKiJiIiIiIj4GCVqIiIiIiIiPkaJmoiIiIiIiI9RoiYiIiIiIuJjlKiJiIiIiIj4GCVqIiIiIiIiPkaJmoiIiIiIiI9RoiYiIiKSTxITEzHGEB8fn23ZjSYuLg5jDMnJyc6y+Ph4jDEkJlPsTj8AACAASURBVCbmW1w5iYqKIiwsLL/DkOuEEjURERERydYDDzyAMYZ27drldygeOXDgAC+88AKdOnUiNDQUYwyxsbFZ1g8LC8MY4/bryJEjLnXj4+OZMmVKHp+BK8uymDt3Lg8++CA1a9YkKCiIatWqce+997Jhwwa3+1y+fJnJkydTp04dAgMDqVq1KsOGDePMmTMu9Y4dO8abb75Jx44dqVq1KkWLFqV27do88cQT/Pe//3Xb9okTJxg8eDCVK1cmMDCQ+vXr8/bbb2NZltfP/Ubmn98BiIiIiMj/REZGkpKSQuHChfM7FACOHj3KkiVLCA8PZ82aNSQnJ+f5qNHo0aN57rnnCAgIyNX+O3fu5NVXX6Vq1ao0adKE5cuX57hPnTp1eOGFFzKVFy9e3OX3+Ph4kpOTefrpp3MVW26cP3+eRx99lIiICB588EFq1KjBgQMHmDlzJi1atODDDz/kkUcecdnnmWeeYerUqXTp0oVhw4aRlJTE1KlT2bZtG6tWrcLPzx6v2bBhA8OGDaNdu3YMGjSIkJAQfvrpJ9555x3mz5/Pd999R7169ZztXrhwgQ4dOrBt2zYGDx5M3bp1Wb58OX//+985ePAgcXFx16xfrndK1ERERER8iJ+fH4GBgfkdhtPcuXO5ePEi8+bNo0WLFsyePZtx48blybFOnTpF8eLF8ff3x98/929TGzVqxKFDhwgNDeXIkSOEhobmuE/58uUzJTu+wt/fn8TERNq0aeNS3q9fP+rXr8+wYcN4+OGHncnX9u3bmTZtGl27dmXhwoXO+jVq1OCpp57i//2//8fDDz8M2Anqzp07CQ8Pd2n7rrvuokOHDowZM4YFCxY4y99//302bdrE1KlTGTx4sDOO+++/n1dffZXHH3+c6tWr50k/3Gg09VFERERuKGnXOq1evZqXXnqJ6tWrU7RoUZo1a8b69esBWLt2La1atSI4OJiKFSvy8ssvu21r8+bNdOnShZCQEAICAqhduzbjx4/n0qVLmep+9tln3Hrrrc5paC+++CIXL17MVM/dNWqXL19m/PjxREZGUqFCBYoUKUK1atUYMGAAR48eddk/OTkZYwxxcXEsXbqUJk2aEBgYSMWKFRkxYoTb2LLzwQcfEBUVRaNGjbj77ruJj4/n8uXLmerFxsZijOHw4cP06tWLsmXLEhwcTLt27di6dWuWMc6bN49GjRpRtGhR5xt/d9eoXYnixYt7lJxldOnSJU6ePJnl9rCwMNauXcuePXtcpkdmvG5u//79PPTQQ5QuXZqgoCBiYmL45ZdfrjieNP7+/pmSNLCTyzZt2nDo0CEOHTrkLP/000+xLCvTqF+/fv0ICgpi7ty5LueUMUkDaN++PWXKlOGnn35yKf/kk08ICgqiX79+LuVPP/20M6EX71CiJiIiIjek5557jiVLljBkyBDGjh3Lrl276NixI0uWLKFr1660bt2a119/nTp16jBmzBiXN7cAX3zxBbfffju//PILw4YNY+rUqbRo0YIxY8bw0EMPudRdvHgxXbp04cSJE4wZM4bBgwezYMECnn32WY9ivXDhApMmTaJWrVqMGDGCqVOn0qFDB2cSdeHChUz7LFu2jN69e3PHHXcwefJkGjZsyOuvv87EiRM97qNNmzbxn//8h8ceewywk7G9e/eyatWqLPfp1KkTBw4cIC4ujqeffprNmzfTpk2bTG/4AZYsWcKAAQPo1KkTU6dO5Y477vA4Nm/bsGEDQUFBlCxZklKlSvHYY4+xf/9+lzpTpkyhTp06hISE8NFHHzm/6tat66xz5swZIiMjKVSoEK+++iqDBg0iMTGR++67j9TUVK/H/ccff1CkSBFKlSrlLNu0aRN+fn40bdrUpW5gYCARERFs2rQpx3ZPnDjBqVOnKF++vLPs8uXLbN261fmBQ3pNmzbFGONR2+IZTX0UERGRG1Jqairr16+nSJEiANSrV4/77ruP7t27s27dOho3bgxAnz59qF69OjNmzHBOjTt37hx9+vShWbNmJCQkOKfpPfnkkzRs2JChQ4eSmJhIVFQUqampDBkyhDJlyrBx40ZCQkKcdW+55RaPYg0ICODAgQMULVrUWda/f39atmxJ3759WbJkCT169HDZZ/v27Wzfvt15PVn//v1p0KAB06ZN4/nnn/fouLNmzSI4OJj7778fgDvuuIPQ0FA++OADOnbs6Haf6tWrs3DhQowxAHTt2pUmTZowfPhwVqxYkSnGH3/80SXRyQ/169enb9++1K1bl4sXL5KYmMj777/P6tWr2bhxI5UqVQKgc+fOTJkyhZSUlCynSR45coQRI0YwcuRIZ1loaCgjR45k1apVxMTEeC3uZcuWsXHjRh599FGXxGn//v3OUd6MKleuzHfffceFCxecj313xo8fz8WLF51JOtgLj6SkpFC5cuVM9QMCAggJCWHfvn1XeVaSRomaiIiIZGnFbyv48/Sf+R2GiwrFKtCpZqerbmfAgAEub1Rbt24NQLNmzZxJGkCRIkVo2rQp3377rbPsq6++4uDBg7z22mscP37cpd0777yToUOHsnLlSqKiotiyZQv//e9/GT58uDNJAyhZsiT9+/f3KGkyxjiTtNTUVE6dOsWlS5eIjo4G7NGgjIla586dXRb9MMbQtm1bpk+fzunTpylWrFi2x0xJSeHTTz/l/vvvd9YtXLgwPXv25O233+avv/6iTJkymfYbOXKkM0kD+3qxDh06sGrVqkzHveuuu/I9SQN7dDS9Bx98kMjISHr27MnYsWN57733PG7Lz8+Pp556yqUs7e/066+/ei1R+/XXX3n00UepXLkyb7zxhsu2s2fPZrkQS1pCd/bs2SwTtQULFvD666/TqVMnHn/8cZd2gWzbTqsjV09TH0VEROSGdNNNN7n8Xrp0acBecCGj0qVLu1wLlpSUBEDv3r0JDQ11+apTpw4ABw8eBGDXrl0AzvL00q+ml5P58+fTrFkzihYtSunSpQkNDXWew7Fjx3I8P4CyZcsCZLquzZ0FCxZw4sQJ2rRpw2+//eb8ioyM5Pz585mmgqZxl3jVq1eP1NRU9uzZ41J+88035xhHfnn44YcJCwvLlMTlpFKlSpmmBV5Jv3ti9+7dtGvXDmMMy5cvz3Q9XlBQEOfPn3e777lz55x13Fm2bBk9e/akUaNGzJs3zyXpTtsnu7azaleunEbUREREJEveGLnyVYUKFbqi8vTS7hc1adIkIiIi3NZJmy7nDYsWLeKBBx6gadOmvPnmm1StWpXAwEBSU1Pp1KmT28U9sjsPT+539cEHHwD21E93Zs2alWnk6Er5+pv6sLAwl5FUT1xtv+ckOTmZtm3bcvr0aVavXk2DBg0y1alUqRI///wz58+fzzT6tW/fPkJCQtyOpq1YsYKuXbtSv359Vq5cSYkSJVy2ly5dmqJFi7qd3nj+/HmOHDnidtETyR0laiIiIiJXqFatWgAEBwfTvn37bOumjWzt2LEj07aff/7Zo+N99NFHBAYGsmbNGpfkxl2b3vD777/z9ddf07NnTzp37pxp++rVq5k5cyZbtmyhUaNGLtuSkpJo3ry5S9nPP/9MoUKFCtyy7b/99pvLYhqAywjTtZacnExUVBQnTpxg1apV3HrrrW7rNWnShJUrV7Jx40bnlF6wR7y+//57IiMjM+2zYsUKOnfuTJ06dVi1apVzhDk9Pz8/brvtNrZt25YpCdy4cSOWZblMG5aro6mPIiIiIlcoJiaGcuXKMWHCBP76669M21NSUjh16hRgX6NVpUoVZs+ezZEjR5x1Tp48ycyZMz06XqFChTDGuIycWZbFK6+8cpVn4t6sWbOwLIuhQ4fSrVu3TF9pq1XOmjUr074TJ050GTnaunUrq1atol27djleF5cf3P39AGbMmMEff/zBPffc41JerFgxjh075pXRsSuxZ88e2rZty/Hjx1m5cmWmBDm9Bx54AGMMU6ZMcSl/7733OHv2LD179nQpX7lyJV26dKF27dqsXr3a7bWHaR566CHOnj3Lu+++61I+ZcoU/P39eeCBB3JxduKORtRERERErlBwcDAffvghnTt3pnbt2vTu3ZuaNWty/PhxduzYwaJFi1i8eDFRUVEUKlSIyZMn06NHD5o2bUq/fv3w9/dn1qxZlC1blr179+Z4vG7durFw4UKio6Pp1asXFy9eZMmSJXmycENqairx8fGEhYVx2223ua0TFhZGo0aN+OSTT3jjjTdcrsnas2cPMTEx3HvvvRw4cIDp06dTtGhRJk2a5PVYs5OWxKb10Y8//ugsi4yMdI4qffjhh3zwwQd06tSJsLAwLl26RGJiIkuWLCE8PDzTzb2bN2/O0qVLGTRoEC1btqRQoUJER0dTrly5K44xLCyMPXv25Jj0nTp1irZt25KcnMzgwYPZuXMnO3fudKnToUMH5+hfgwYNGDhwINOnT6dr167ceeedJCUlMXXqVNq0aeO82TXY9wK87777sCyLxx9/nOXLl2c6fvoVLvv168fs2bMZOnQoycnJ1K1bl2XLlrF48WJGjx7tsoCNXB0laiIiIiK5EBMTw6ZNm5gwYQJz587l8OHDlC5dmvDwcIYOHeqy9H63bt1YsGABL730EnFxcZQrV47Y2FgiIyOzXOY+vQcffJBTp04xefJkhg8fTunSpbnnnnuYMGGCc6EKb1mxYgX79+9n6NCh2da7//77ef7551m0aJHLG/8VK1YwdOhQxo4dS0pKCs2bN2fSpEke34rAW1588UWX37dt28a2bdsAGDt2rDNRa9KkCQkJCcybN4/Dhw9jWRY1atTg2Wef5bnnnnO5PxnAM888w65du1iwYAEzZ87k8uXLrFmzJleJ2unTpz26lvHo0aPs3r0bgGnTprmts2bNGpdpmlOmTCEsLIx3332XL774gpCQEAYPHsxLL72En9//JtX99NNPzgVGnnnmGbdtp0/UihQpwqpVqxg9ejSffvopR48eJTw8nGnTpjFw4MCcT1o8Zq5k2LZx48bW5s2b8zAc73l/3YKr2r9vi24e142KjwIgMTbxqo4ptvjv4wGIjYjN1zhERG4kSUlJPrFMuhRcsbGxzJkz55pPCSyofvzxRxo2bMisWbNclsCX/OXJc6ExZotlWXl+MZ6uURMRERERuca+/PJLGjZs6HJDaZH0lKiJiIiIiFxjI0aM4Pvvv3eZhiiSnh4ZIiIiIiIiPkaJmoiIiIhctfj4eF2fJuJFStRERERERER8jBI1ERERERERH6NETURERERExMcoURMREREREfExStRERERERER8jBI1ERERERERH6NETURERERExMcoURMREREREfEx/vkdQFbeX7cgv0MQERER8ZrExETatm3L7NmziY2Nze9wciUuLo5x48axe/duwsLCcqxvjOGxxx4jPj4+z2PLT9fD31Z8j88maiIiIuIbfO3D074tuuV3CCJ5YtmyZbzyyiv88MMPBAQE0K5dOyZOnEiNGjXyOzTJB5r6KCIiInINREZGkpKSwqOPPprfoeTa6NGjSUlJoXr16vkdynVn0aJF3H333aSkpDBp0iRGjBjB119/ze23387+/fvzOzzJBxpRExERkRteamoq58+fJygoKM+O4efnR2BgYJ61fy34+/vj76+3j9528eJFBg8eTNWqVfnmm28oVqwYAHfccQeNGjUiLi6Od999N5+jlGtNI2oiIiJyQ4mPj8cYw6pVq3j55ZcJDw8nMDCQ+fPnO+ts3ryZLl26EBISQkBAALVr12b8+PFcunQpU3sLFy6kYcOGBAYGUq1aNcaNG8eqVaswxrhcm5WYmJipDODMmTOMGjWK8PBwAgICqFChAr169WLPnj0u9dLvP3v2bOrXr09AQADVq1dn4sSJHp37nj17MMYwduxYl/KYmBiMMUyePNmlvFmzZtStW9f5e1xcHMYYkpOTXept376dTp06ERwcTJkyZejZsyeHDh3KMo558+bRqlUrihcvTlBQEM2aNWPBgpyn2J4/f56iRYvy2GOPuZQ/+eSTGGMYMmSIS/kDDzxAiRIlXP5uJ06c4Nlnn6VmzZoEBAQQGhrKQw89xK5du9we79VXX6V+/foEBgZSqlQp7rnnHrZt25ZjrABz5syhcOHCdOvWjXPnzmVZb+3atezfv5++ffs6kzSAiIgIoqKimDdvHhcvXvTomHL90EciWbiS+fgHTh52u4/m0IuIiPiu4cOHc/HiRfr160eJEiWoXbs2AF988QVdu3alZs2aDBs2jDJlyrBu3TrGjBnD999/z//93/8525g3bx4PPfQQ4eHhjB07Fn9/f+bMmcPnn3/uUQwXL14kJiaGb7/9lm7dujFs2DB+/fVX3n77bVauXMnmzZupUqWKyz4zZ87k4MGD9OnTh1KlSjF37lyeffZZqlSpwsMPP5zt8apXr85NN91EQkIC48aNA+DChQv8+9//xs/Pj4SEBJ555hkATp48yZYtW3jyySezbXP37t20bt2a8+fPM2jQIKpWrcrnn39Op06d3NYfPXo048ePp1OnTrz88sv4+fmxePFiunfvzvTp0xk4cGCWxwoICKBly5asWbPGpXz16tXO+NNYlkViYiKtW7d2jgKeOHGCli1bsnfvXnr37k39+vU5cOAAb731Fs2aNWPz5s3OaZ0XL16kU6dOfPfddzz66KMMGjSIEydO8N5773H77bfz9ddf07hx4yxjffXVV3nhhRcYOHAgU6dOxc8v6/GRTZs2AdCiRYtM25o3b05CQgK//PIL9evXz7INuf4oURMREZEbUkpKCtu2bXOZ7nju3Dn69OlDs2bNSEhIcL7Bf/LJJ2nYsCFDhw4lMTGRqKgoLl26xNChQwkNDWXjxo2ULl0agAEDBnDLLbd4FEN8fDzffvstI0aMcBkVa9++PXfffTejRo3io48+ctln7969JCUlUbJkSQB69+5N9erVmTZtWo6JGkB0dDRz5szh7NmzBAUFsX79es6ePcsjjzzCZ599xqVLl/D392ft2rWkpqYSHR2dbXsvvPACx44dIyEhgbZt2wIwcOBAunbtmmnkaevWrYwfP55Ro0bx6quvOsufeuopOnfuzKhRo+jVqxfFixfPNv6EhAR+/fVXatWqxd69e/n999955JFHmDt3LgcPHqR8+fL89NNPHDp0yCX+MWPGsGvXLtavX0/Dhg2d5bGxsTRo0ICxY8c6RzynT59OYmIiK1asICYmxln373//O3/7298YPnw4iYmJmeK7fPkygwcP5q233mL8+PE8//zz2fYf4LwGrXLlypm2pZXt27dPidoNRlMfRURE5IY0YMCATNekffXVVxw8eJDHH3+c48ePc+TIEefXnXfeCcDKlSsB2LJlC/v37yc2NtaZpAEUK1aM/v37exTD4sWL8fPzY9SoUS7ld911FxEREXz22WdcvnzZZdvjjz/uTNIAgoKCaN68Ob/++qtHx4yOjubixYt88803ACQkJFCuXDmGDBnCqVOnnKM7a9aswRjjTL7cuXz5Mp9//jmNGzd2qWeMYeTIkZnqf/zxx84l+9P37ZEjR7j33ns5deoU69atyzH+tLjTvhcqVMg5LTOtPG3ULa2+ZVl8/PHHREZGUrlyZZdjBwcH07x5c+ffFmDu3LnUqVOHRo0audS9cOECHTp04N///jcpKSkusZ07d45u3brx7rvvEh8f71GSBnD27FnAHjHMKO26xrQ6cuPQiJqIiIjckG6++eZMZUlJSYA9SpWVgwcPAvaUP8A5ZTI9d2Xu7N69m0qVKrkkemnq16/P999/z5EjRyhXrpyz/KabbspUt2zZshw9etT5+4kTJzIlEaGhoRQqVMgl0YmJiXGOhN12222ULl2ahIQEWrRoQUJCAg0bNqRMmTJZxn/o0CFOnz5NnTp1Mm2rV69eprKkpCQsy3JbP01a/2alSZMmFC9enISEBJ588kkSEhJo3Lgx4eHhNGjQgISEBB566CESEhIoU6YMERERABw+fJijR4+ycuVKQkND3badfnpiUlISKSkpWdYFOHLkCFWrVnX+PnLkSE6dOsXHH3/s0ehmmrQPDM6fP59pW9q1bXm50I34JiVqIiIickNy98bXsiwAJk2a5HyDn1GlSpXyNK6cFCpUKMc6Q4YMYc6cOS5laTepLl++PPXq1SMhIYGzZ8+yYcMGpk2bhp+fH23atGH16tX079+fH3/80Xm9mrdYloUxhuXLl2d5HjlN7/P396d169asWbMGy7JISEigV69egD16ljYKuXbtWqKjozHGOI8N9rTSZ5991qNYGzRowD//+c8s62RM4jp37szChQuZNGkSMTExlC1bNsfjwP8eU/v27XNZvCWtDNxPi5TrmxK1PHS1NwjVYiQiIiLXVq1atQAIDg6mffv22dYNCwsDYOfOnZm2uStz56abbmLFihUcP36cUqVKuWz7+eefKVGiBCEhIR61ld7IkSN55JFHXMoqVKjg/Dk6Opq33nqLzz//nAsXLtCuXTsA2rVrx/Dhw1m+fDmWZeV4fVpoaCjFihVjx44dmbb9/PPPmcpq1arFihUrqFatWqaE5EpER0ezbNkyFixYwL59+1zinzJlCosWLeL48eMu8YeGhlKqVClOnjyZ4982LdbDhw8THR2d7UIgGePq3bs3d999N23btmXVqlUuo6FZadKkCQDr1q3LFNv69espUaKE2xFgub7pGjURERERh5iYGMqVK8eECRP466+/Mm1PSUnh1KlTADRu3JiKFSsSHx/PsWPHnHVOnz7NzJkzPTpe586duXz5MhMmTHApX758Odu2bePee+/1OElIr169erRv397lK/093KKjo7l8+TLjxo2jWrVqhIeHO8vPnz/Pa6+9hr+/P5GRkdkep1ChQtx9991s3rzZZSVGy7Lc3jIg7Wbfzz//PKmpqZm25zTtMX38AGPHjiUgIIDbb78dsG8qXqhQIeftB9Inan5+fvTs2ZONGzdmeSuA9LcU6NWrF3/++WeWI2pZxRoVFcWKFStITk6mbdu2/PnnnzmeT5s2bahYsSLvv/8+p0+fdpb/8MMPJCYm0r17dwoXLpxjO3J90YiaiIiIiENwcDAffvghnTt3pnbt2vTu3ZuaNWty/PhxduzYwaJFi1i8eDFRUVH4+/vz+uuv07NnT5o2bUqfPn3w9/cnPj6esmXLsnv3bue0u6zExsYyZ84c/vGPf5CcnExkZCS//fYbb731FuXLl3dZGdGboqKi8PPzIykpidjYWGd5vXr1qFChAj///DPNmzfPdvXFNK+88grLly/n7rvvZvDgwVSpUoXPP/+cw4cPZ6rbpEkT4uLiiIuLIyIigu7du1OpUiUOHDjAli1bWLZsGRcuXMjxmBEREZQpU4akpCSioqKcSWiJEiVo3LgxGzZsoGLFiplG7caPH8+3335Ljx496NGjB82bN6dIkSLs2bOHZcuW0ahRI+eqj0OGDOGrr75ixIgRJCQkEB0dTYkSJdi7dy+rV68mMDAw020C0rRq1YqVK1fSqVMnoqKiSEhIyHbKbOHChXnzzTd54IEHaN26Nf369ePkyZNMnjyZ0NBQ560U5MaiRE1EREQknZiYGDZt2sSECROYO3cuhw8fpnTp0oSHhzN06FCXpfcffvhhChcuzMsvv8zYsWMpX748ffr04ZZbbqFr164ULVo022MVLlyYL7/8kldeeYV58+axaNEiSpUqRffu3XnllVdcFqrwptKlSxMREcHWrVszTW+Mjo7mk08+yXHaY5rw8HC++eYbhg0bxrRp0wgICOCOO+7go48+onz58pnqjx07lsaNGzN16lSmTJnCmTNnKFeuHH/729+YOnWqR8c0xhAVFcWiRYsyxdmuXTs2bNjgdrXKkiVL8u233/LGG28wf/58PvvsM/z9/alSpQqtWrWib9++zrqFCxfmiy++4K233uKjjz5yjtJVqlSJpk2bZrrpdkbNmzdn1apVdOzYkTZt2pCQkJDt37N79+4ULVqUV155heHDhxMQEEC7du34xz/+oevTblAm7cJKTzRu3NjavHlzHobzP1d7fde1NGmT/Y87ool3P+24Ua9Ri/8+HoDYiNh8jUNE5EaSlJR0VdcMias33niD4cOHs27dOpo3b57f4YiIhzx5LjTGbLEsK+u7nXuJrlETERERyaULFy5kutbq9OnTzJgxg7Jly3LbbbflU2QiUtBp6qOIiIhILu3atYs77riDBx98kBo1anDgwAHmzJnD7t27efvttylSpEh+hygiBZQSNREREZFcCg0NpXnz5nz88cccOnQIf39/GjRowIQJE+jRo0d+hyciBZgSNREREZFcKlu2LJ9++ml+hyEi1yFdoyYiIiIiIuJjlKiJiIiIiIj4GCVqIiIiIiIiPkaJmoiIiIiIiI9RoiYiIiIiIuJjlKiJiIiIiIj4GCVqIiIiIiIiPkaJmoiIiIiIiI9RoiYiIiJyDSQmJmKMIT4+Pr9DybW4uDiMMSQnJ3tU3xhDbGxsnsZ0rUVFRREWFuZSFhsbizEmfwLy0PX4t7je+ed3ACIiIuLbfjz6Q36H4OKWsg3zOwTxEampqVSrVo39+/fz0ksv8eKLL+Z3SDnauHEjc+fOZcuWLfzwww+cOXOG2bNnu02ikpOTqVGjhtt26tevz08//eRSFhcXR0REBJ07d86L0N06duwYH374IV988QVJSUkcOXKEatWq0aZNG1588UWqVq2aaZ8TJ04wevRoFi1axNGjRwkPD2fQoEH079/fJeH95ZdfmDt3LitXruT333/n3LlzhIeH0717d55++mmCg4Mztb1z506effZZ1q5dy4ULF7jtttsYN24c0dHRedoPeUGJmoiIiMg1EBkZSUpKCoULF87vUHJt9OjRPPfccwQEBOR3KAAsX76c/fv3Ex4eTnx8PKNHj87zka2VK1diDfcRiAAAIABJREFUWVau91+2bBkzZsygTp06NGzYkO+++y7Hfbp06ULXrl1dykqVKpWp3rhx43jssceuaaK2YcMGhg0bRrt27Rg0aBAhISH89NNPvPPOO8yfP5/vvvuOevXqOetfuHCBDh06sG3bNgYPHkzdunVZvnw5f//73zl48CBxcXHOurNmzWLGjBnce++99OzZk8KFC7NmzRpGjx7N/PnzWb9+PUWLFnXW//3332nZsiX+/v6MHDmSkiVL8t577xETE8Py5ctp3779NesXb1CiJiIiIje81NRUzp8/T1BQUJ4dw8/Pj8DAwDxr/1rw9/fH39933j5+8MEHhIeH889//pP77ruPxMRE2rZt6/XjpH98FClS5KraGjBgACNGjCA4OJgFCxZ4lKjdcsstPPLII1d13LxSp04ddu7cSXh4uEv5XXfdRYcOHRgzZgwLFixwlr///vts2rSJqVOnMnjwYAD69evH/fffz6uvvsrjjz9O9erVAejWrRujRo2iZMmSzv379+9PrVq1GD9+PB988AGDBg1ybhs1ahTHjx9ny5YtREREANCrVy/q16/PwIED2bFjh89PUU1P16iJiIjIDSU+Ph5jDKtWreLll18mPDycwMBA5s+f76yzefNmunTpQkhICAEBAdSuXZvx48dz6dKlTO0tXLiQhg0bEhgYSLVq1Rg3bhyrVq3KdD1aVteonTlzhlGjRhEeHk5AQAAVKlSgV69e7Nmzx6Ve+v1nz55N/fr1CQgIoHr16kycONGjc9+zZw/GGMaOHetSHhMTgzGGyZMnu5Q3a9aMunXrOn/P6hq17du306lTJ4KDgylTpgw9e/bk0KFDWcYxb948WrVqRfHixQkKCqJZs2Yub+Y9cfDgQZYuXUqvXr24887/396dx0dV3f8ff38g7IiyK8oioFApq4jQqmBEqYpI+UkLAgICdasVEa0oyqqiRaVQqD8rCKIo2EItFVAB6RcrKpvypYJWlipqlUWEsCXA+f5x78SZ5E4mIRPmJryej8c8AmfOPefMJyeZ+eTce+41qlWrlqZPnx5Yt0GDBurUqZPWrVun9PR0Va5cWdWqVVP//v1zjTPR/Ai6Rq0gateuHXjKXiKHDx/WwYMHA5/bvn17dgIya9YsmVn2I6dVq1apY8eOqlSpkqpXr67BgwcrIyOjwOOJaNCgQa4kTZI6d+6satWq5To9c86cOapYsaKGDBkSUz506FBlZWVp7ty52WVt27aNSdIifvnLX0pSTNsHDhzQ3/72N3Xq1Ck7SZOkypUra/Dgwfr000+1evXqE3uRKUKiBgAATknDhw/XK6+8oiFDhuj3v/+9mjRpIkl6/fXX9dOf/lSffvqp7rnnHk2ePFkdOnTQww8/rN69e8e0MXfuXPXs2VMHDx7UqFGjdOedd+rVV1/V/fffn68xZGVlqUuXLpowYYLatGmjp59+Wr1799a8efN08cUXa8eOHbmOeeaZZzR27Fj17t1bTz75pM466yz99re/1Zw5cxL2V79+fTVs2FDLly/PLsvMzNQ777yjUqVKxZTv27dPa9euTXhtz7Zt23TppZdq5cqV+vWvf62xY8dq165d+tnPfhZYf+TIkerVq5dOO+00jRs3ThMmTFDFihXVs2dPTZ06NeFriHjhhRd07Ngx3XTTTUpLS1OfPn00f/58ff/994H1d+zYoSuuuEINGzbUE088oR49emj27Nm6/PLLAxOgePMjFZ588klVrFhRlSpVUt26dfXwww/ryJEj2c/XrFlTs2fPliRdeumlmj17dvYj2ocffqiuXbvqoosu0lNPPaWrrrpK06dP17Bhw5I+5u+//1779+9X7dq1s8uOHz+udevWqXXr1rlWl9u1ayczy1cyFfm5iG57w4YNOnLkiDp06JCrfvv27SWp2CVq4Vm7BgAAOIkOHTqk9evXx5zuePjwYQ0aNEgXX3yxli9fnn2a3y233KKWLVtq2LBhWrFihTp16qSjR49q2LBhqlmzpj744ANVrVpVkndqW4sWLfI1hpkzZ+qf//yn7r333phVsc6dO6tr164aMWJErg/bn3/+uTZt2pS90nDzzTerfv36mjJlim688caEfaanp2vWrFk6ePCgKlasqPfee08HDx5U37599dprr+no0aNKS0vTP/7xDx07dixhovbggw/qu+++0/Lly7NPO7zjjjvUo0cPrV+/PqbuunXr9Mgjj2jEiBF69NFHs8t/85vfqHv37hoxYoRuuukmnXbaaQlfx4wZM3TZZZdlr271799fTz/9tObMmaPbbrstV/0tW7bo6aef1tChQ7PLmjVrpmHDhmny5Mm5kuug+XGylSpVSunp6erevbvq16+vnTt3at68eRo3bpxWrVqlJUuWqHTp0qpUqZL69u2rfv36qWHDhnFPk9ywYYNWrVqliy++WJI3r/ft26fnn39eTz31lCpXrpy0sT/yyCPKyspS//79s8u+++47HTp0SGeffXau+uXKlVONGjX05Zdf5tnusWPHNG7cOKWlpcXM96+++kqSAtuOlCVqO2xYUQMAAKek2267LdeH8LfeekvffPONBg4cqL1792rXrl3Zj2uuuUaSt5mEJK1du1ZfffWVBgwYkJ2kSd6pVrfeemu+xrBgwQKVKlVKI0aMiCm/9tpr1apVK7322ms6fvx4zHMDBw6MOR2sYsWKat++vf7973/nq8/09HRlZWVp5cqVkqTly5erVq1auuuuu7R///7sVYe3335bZpbnNV/Hjx/XwoUL1bZt25h6Zqb77rsvV/2XXnpJZqb+/fvHxHbXrl3q1q2b9u/fr1WrViV8De+++642b94ckwS0bNlSrVq10owZMwKPqVKlim6//faYsttvv11VqlTRggULctUPmh8nW7169bRs2TLdeeed6tatmwYNGqQ33nhDQ4YM0dKlS/XKK68UqL0OHTpkJ2kR6enpOnr0aL5vuZAff/7znzVx4kT97Gc/08CBA7PLIyuX8TajKV++fNzTOyOGDh2qVatWaezYsTGrnHm1HVm9S9R22JCoAQCAU9L555+fq2zTpk2SvFWqmjVrxjyaNm0qybs2SvJO+ZMUeEpcfk+T27Ztm+rUqROT6EU0a9ZM+/fv165du2LKGzZsmKtu9erVtXv37uz/f//99/rvf/8b8zh27JgkZa+QRU5zjKyEtWnTRlWrVo0pb9mypapVqxZ3/N9++60yMjKyYxMteqe/iE2bNsk5p6ZNm+aK76BBgyT9EN+8TJ8+XWXKlFHr1q312WefZT+6dOmiNWvWaMOGDbmOadiwYa6NQMqVK6eGDRtq69atueoHzY+wePDBByV5p+kWRLy5Iylm/hTGokWL1KdPH1144YWaO3duzHVykcQ3+rTNaIcPH84zOX7ooYf0hz/8Qb/61a9y/XEjr7YPHz4cU6e44NRHAABwSgr60BbZdv13v/tdzIYE0erUqVOk40qkdOnSCevcddddmjVrVkzZtm3b1KBBA9WuXVsXXHCBli9froMHD+r999/XlClTVKpUKXXs2FHLli3Trbfeqg0bNujuu+9O6tidczIzLV68OO7raNasWZ5tZGRkaN68ecrKylLr1q0D68yYMUOTJk0q1FjD/KG+bt26Kl26dK4kPpG85k5hbjkQsWTJEvXo0UPNmjXTm2++qSpVqsQ8X7VqVVWoUCHwFMQjR45o165d6tixY2Dbo0eP1vjx4zVw4EA988wzuZ6P/FwGtR0pCzotMsxI1AAAAHznnXeeJKlSpUoJ77kUuTbqk08+yfVcUFmQhg0basmSJdq7d2+u+2J9/PHHqlKlimrUqJGvtqLdd999ua5TOvPMM7P/nZ6ermnTpmnhwoXKzMzUFVdcIUm64oorNHz4cC1evFjOuYTXp9WsWVOVK1fW5s2bcz338ccf5yo777zztGTJEtWrVy9mN8mCmDdvnjIyMvToo49mf7+iTZ48WS+++KKeeOKJmBW0rVu3KjMzM6bsyJEj2rp1a+CKYJht3bpVx44di9lMI9WWLFmi7t27q2nTplq6dGngKnGpUqXUpk0brV+/XkeOHIk5TfGDDz6Qc05t27bNddzo0aOz7xH33HPPBe5m2bx5c5UrVy7w1Nn33ntPkgLbDjNOfQQAAPB16dJFtWrV0oQJE7Rnz55czx86dEj79++X5H3oO+usszRz5kx999132XUyMjIC/+IfpHv37jp+/LgmTJgQU7548WKtX79e3bp1U6lSBf+4dsEFF6hz584xj+hd9tLT03X8+HGNGTNG9erVy95ePT09XUeOHNFjjz2mtLQ0XXbZZXn2U7p0aXXt2lVr1qzR22+/nV3unAu8ZUC/fv0kSQ888ED2qZjR8nvaY7Vq1XTvvffqhhtuyPUYNGiQdu/erddeey3muH379mnatGkxZdOmTdO+fftO6g2iCyLodMTjx49r5MiRkqTrrrsu5rnKlSsHztui9uabb+rnP/+5mjRpomXLluV5umzv3r118OBBPfvsszHlkyZNUlpaWvbW+xFjx47VmDFj1K9fP82YMSPuz0PlypV13XXXacWKFfroo4+yyzMyMvTcc8/pvPPOU7t27QrxKk8+VtQAAAB8lSpV0gsvvKDu3burSZMmuvnmm9W4cWPt3btXmzdv1vz587VgwQJ16tRJaWlpmjhxovr06aN27dpp0KBBSktL08yZM1W9enVt27Yt4c11BwwYoFmzZunxxx/X9u3bddlll+mzzz7TtGnTVLt27ZidEZOpU6dOKlWqlDZt2qQBAwZkl19wwQU688wz9fHHH6t9+/b52n1x/PjxWrx4sbp27ao777xT55xzjhYuXKidO3fmqnvRRRdp9OjRGj16tFq1aqWePXuqTp06+vrrr7V27VotWrRImZmZcfvavHmz3n33XQ0YMCDujbe7deumMmXKaPr06erZs2d2eaNGjTRmzBht3LhRF154odauXasZM2aoadOm+s1vfpPwdSbLf/7zn+ydPP/1r39JkhYuXJi95Xy/fv2yb/g8ZMgQ7du3Tz/5yU9Ut25d7dq1S3/5y1+0du1aXX/99brhhhti2m7fvr2WLl2qxx9/XPXq1ZOZqVevXgUe4/bt23XuueeqY8eOWrFiRZ5116xZo+uvv17OOQ0cOFCLFy/OVSd6dXfIkCF6/vnnNWzYMG3fvl0/+tGPtGjRIi1YsEAjR46MuUfd1KlTNWrUKNWrV0+dO3fOdQuK2rVr68orr8z+/2OPPaZly5bpqquu0t13360qVaroT3/6k7788ku9/vrrxepm1xKJGgAAQIwuXbpo9erVmjBhgl588UXt3LlTVatWVaNGjTRs2LCYrfdvvPFGlSlTRuPGjdOoUaNUu3ZtDRo0SC1atFCPHj1UoUKFPPsqU6aM3njjDY0fP15z587V/PnzdcYZZ6hnz54aP3686tatWySvsWrVqmrVqlX2DaCjpaena86cOQlPe4xo1KiRVq5cqXvuuUdTpkxRuXLldPXVV2v27NmBp+aNGjVKbdu21eTJkzVp0iQdOHBAtWrV0o9//GNNnjw5z74iN7Tu0aNHnq/t8ssv11tvvaUvvvgiO4bnnHOO5s2bp+HDh+vll19W2bJl1adPH02cOPGEbkB9orZt26aHHnoopmz+/PmaP3++JOmSSy7JTtSuvfZazZ49W88++6z27NmjcuXKqVmzZpo6dapuvfXWXKtL06ZN0x133KFHHnkke+X3RBK1yLH5uaZr48aN2Zt1xLumMTpRK1u2rJYuXaqRI0fq5Zdf1u7du9WoUSNNmTJFd9xxR8xxkR1IP//885gdPiM6duwYk6g1btxY//znP3X//fdrwoQJyszMVJs2bbRkyZKEpzKHkRXkwsG2bdu6NWvWFOFwfvDcqoLdnT6Vfrd6lCTp3ovGJLXdwR1uSFypBJr54UxJ0oBWA1I6DgA4lWzatOmErxlCbk8++aSGDx+uVatWZd9sF6nToEEDNWjQIOHqEDyTJ0/W8OHDtXHjxlDvflkU8vO70MzWOueK/II3rlEDAAA4QZmZmbmutcrIyNDUqVNVvXp1tWnTJkUjA07cG2+8oVtuueWUS9LChlMfAQAATtDWrVt19dVXq1evXjr33HP19ddfa9asWdq2bZv++Mc/5rpvF1AcFPT+bCgaJGoAAAAnqGbNmmrfvr1eeuklffvtt0pLS1Pz5s01YcIE/eIXv0j18AAUYyRqAAAAJ6h69ep6+eWXUz0MJLB9+/ZUDwEoMK5RAwAAAICQIVEDAAAAgJDh1McQK+wtCk7V7f0BAACA4o4VNQAAIEkqyL1VAaCkCdvvQBI1AACgsmXL6tChQ6keBgCkzKFDh1SuXLlUDyMbiRoAAFCNGjW0Y8cO7dmzR1lZWaH7yzIAFAXnnLKysrRnzx7t2LFD1atXT/WQsnGNGgAA0Omnn65y5cpp586d2r17t44ePZrqIQHASZGWlqby5curXr16Kl++fKqHk41EDQAASJLKly+vunXrpnoYAABx6iMAAAAAhA6JGgAAAACEDIkaAAAAAIQMiRoAAAAAhAybiSCUvsj4XBt2f5TqYcTVonrLVA8BAAAAJRgragAAAAAQMiRqAAAAABAyJGoAAAAAEDIkagAAAAAQMiRqAAAAABAyJGoAAAAAEDIkagAAAAAQMiRqAAAAABAyJGoAAAAAEDIkagAAAAAQMiRqAAAAABAyaakeAMLruVV/LnQbgzvckISRAAAAAKcWVtQAAAAAIGRI1AAAAAAgZEjUAAAAACBkSNQAAAAAIGRI1AAAAAAgZEjUAAAAACBkSNQAAAAAIGS4j1oRyyy9P2V9b9j9UZ7Pt6je8iSNBAAAAEBBsKIGAAAAACFDogYAAAAAIUOiBgAAAAAhQ6IGAAAAACFDogYAAAAAIUOiBgAAAAAhQ6IGAAAAACFDogYAAAAAIcMNr09hiW6InYybdSfqI8gXGZ8Xul8AAACgOGNFDQAAAABChkQNAAAAAEKGRA0AAAAAQoZEDQAAAABChkQNAAAAAEKGRA0AAAAAQoZEDQAAAABChkQNAAAAAEKGRA0AAAAAQoZEDQAAAABChkQNAAAAAEKGRA0AAAAAQoZEDQAAAABChkQNAAAAAEImLdUDQMn2zpZ1BT7m02//4x1bdp0uadQm2UMCAAAAQo8VNQAAAAAIGRI1AAAAAAgZEjUAAAAACBkSNQAAAAAIGRI1AAAAAAgZEjUAAAAACBkSNQAAAAAImRJxH7XM0vtT2r+zY6EYBwAAAICSgRU1AAAAAAgZEjUAAAAACJkSceojgr2zZV2qhwAAAADgBLCiBgAAAAAhQ6IGAAAAACFToFMfDx09pA27PyqqscRgB0XgxJysn9ET0aJ6y1QPAQAAoFhgRQ0AAAAAQoZEDQAAAABChl0fgRMQ5tMLAQAAUPyxogYAAAAAIcOKGoCTJswrkWx0AgAAwoQVNQAAAAAIGRI1AAAAAAgZTn1EqL2zZV2hjr+kUZskjQQAAAA4eVhRAwAAAICQIVEDAAAAgJAhUQMAAACAkCFRAwAAAICQIVEDAAAAgJAhUQMAAACAkCFRAwAAAICQIVEDAAAAgJAhUQMAAACAkCFRAwAAAICQSUv1AICS7p0t6wp1/CWN2iRpJAAAACguWFEDAAAAgJAhUQMAAACAkCFRAwAAAICQ4Ro1AJC0YfdHqR5CsdSiestUDwEAgBKJRA1IoLCbgQAAAAAFRaIGACiRwrxKykokACARrlEDAAAAgJAhUQMAAACAkOHURwDACQvz6YUAABRnrKgBAAAAQMiYcy7/lc12SvpP0Q0nKWpI2pXqQeRDcRlncUAsk4t4JhfxTC7imVzEM3mIZXIRz+QinsnVxDl3WlF3UqBTH51zNYtqIMliZmucc21TPY5Eiss4iwNimVzEM7mIZ3IRz+QinslDLJOLeCYX8UwuM1tzMvrh1EcAAAAACBkSNQAAAAAImZKYqD2b6gHkU3EZZ3FALJOLeCYX8Uwu4plcxDN5iGVyEc/kIp7JdVLiWaDNRAAAAAAARa8krqgBAAAAQLFGogYAAAAAIROaRM3MSpnZ3Wa22cwOm9kXZvakmVVK9vFmtsLMXJxHnluXJmGcI8zsVTPb6ve3PUH9i81sqZntN7N9ZrbEzFrlp6/ioDDxNLPzzWysmb1nZjv9GH1oZg/GO97MmpjZX83sOzM7YGYrzSw9+a/s5CtkLJuY2UtmtsnMvjezg347T5nZWXkcUyJjKRX+Zz1HWxWjfub/EKcO8cz7+Hi/szPi1CeeiduoZmYTzewzv42dZva2mV0aUJf3ovjHjs5jfjozywo4psTOzyT8rFc2swfM7H/9+bbLzN41swFmZgH1mZt5H1/bzJ7xj8s0s8/N7Pdmdkac+iV2bkoF/xyeRzs3mdl6MztkZt+Y2XNmFngbs0LNUedcKB6Sfi/JSZovaYikpyRlSVouqVQyj5e0QtJOSX0DHtWKeJxO0m5Jb0naI2l7HnXbSzosaYuku/3HFkn7JTVP9fcs1d93SRP8WLwk6U5Jt0qa67f3kaQKOeo38mP/jaQRkm6XtN7vr3OqY5HiWF7h13vUj8uvJE2RlCHpK0m1TqVYFjaeAW1N9Oeqk/SHgOeJZ+LjnaT/Cfid/UvieULxrC9pm7z3wgmSbvbfY56X1CtHXd6L8j62RcC87CvpiUibp9L8LGQsS0laKemYpBny3ouGSnrfb/Nx5maB4llL0nZJmfLe02/xv2b6c67iqTQ3/deY78/hebRxt9/OCn+OjpX3eelfkiolc46mPGD+i2gm6bikv+Qov9MPxI3JPN4P7Il8Ywo1Tr9uw6h/b8xrHJI+kLRP0tlRZWf7ZW+m+vsWgu97W0mnB5SP94//dY7yef4v/1ZRZZUl/UfSJ/I31ymOj2TMzTjt9vSPv+9UiWWy4ympjaSjkoYpfqJGPBO34STNzGd/xDNxGyslfSHprHzU5b3oxNr9//7x1+YoL7Hzs7CxlNTBr/d0jvKykrZK2pujnLmZ9/GT/Hq9c5T39stH5igvsXMz6vXk+3N4nONrSDrgz73SUeXX+TF9IEf9Qs3RlAfMH3Dkg/WlOcrL+8FYlMzj5Sdq8v5yUyW/E6+w4wxoL+4EkdTY72t6wHPT/R/cM1P9vUvl9z2Pdpv77T4TVVZJ3l80lgXUf8iv3y7VMQlhLNv57T52qsQymfGUVFrSWkl/l9RAAYka8cxfPP3jZ8r7wFY5j3rEM/F75mX+8Xf6/y+jHH9Zj6rLe9EJ/P705+H38pLh0jnKS+z8TMLc7OIff2/Acx9I+jLq/8zNxPH8SNJB5ficK+/z7yFJW6LKSvTcjBOfE0nUBvux6Bfw3BZJH0f9v9BzNCzXqF0kb7AfRBc65w5L+tB/PtnHny1vmfJ7SRlmNt/MmhbxOAsi0taqgOfek2SSLkxif6lQVPE8x//6TVRZC0nlFD+ekfEUV0mJpZmVN7MaZnaOmV0l7y/CkrQoqlpJj6WUvLl5t6Smkn6dRx3imf/Xd4O8Dx37zexbM5tiZqfnqEM8E7++a/yvn5vZQnkf2A6Y2adm1jegL4n3ooLqKe8PwTOdc8eiykv6/CxsLD+QtFfSfWbW08zqmVlTM3tM3jwbnaMvibmZl3KSDjs/M4g6/ri8n/uGZlbDLy7pczNZEs27pmZWOZ91E87RsCRqdSTtcs4dCXjuS0k1zKxsEo/fJu/c8YHyfplOk3S1pPfNrHkRjrMg6kS1G9SX5CWbxVnS42lmpeX95eeopDk5+oq0G9SXVLzjmaxYDpZ3zcoXkt6QdIakvs65lTn6irQb1JdUvGMpJSGeZnaupDGSxjrntifoK9JuUF8S8ZS8Dyqj5SVr/eVdn/FrSSuj3hQjfUXaDepLIp5N/K9/klRNXjxvlnfdymwzG5ijr0i7QX1JxDPIIHl/SZ8R0Fek3aC+pOIdz0LF0jn3naRu8q4dmifvlLtNku6Q9P+cc3/K0Vek3aC+pOIdS6nwc/Nfkqrm3LjC/39V/7/1ovqKtBvUl1T845kMieJkUXUKHdO0go6uiFSUFDQJJW8ZNlInMxnHO+cG5qjzZzP7m7xTIp+SdGURjbMgKvpfg/o7nKNOcVUU8Zwk7xz3B5xzn+ToS3H6KwnxTFYs/ypps7xz0lvLe8OskaNOSY+llJx4PiPvmoqn8tGX4vRHPH3OuYtzFL1gZhskPSLpLv9rpB3F6Y94ek7zv+6XdLlzLlOSzOyv8ubso2Y2y/+rO/H8oU6+3ovMrImkS+SdQrYtoC/F6a8kxDMZscyQd0ra3yS9K++PCXdImmNm1zvn3opqR3H6KwmxlAofz0mSukuaZ2ZD5cW1mV+eJf+056h2FKe/khLPZChInAod07CsqB2Ut9wapHxUnaI6Xv6Kwf9IutzMKhRVPwUQaSeov2T3lSpJjaeZjZP3F/ZnnXOPBfSlOP2VhHgmJZbOuR3OuaXOub8650bJ+0v7E2Y2IkdfitNfSYilVMh4+qePXSnpNudcrq25A/pSnP6IZ95+J+8DyrU5+lKc/oin55D/9eVIkiZlr2b8TdKZ+mHVjXj+UCe/Bvlfn4vTl+L0VxLiWdjfnc3lJWdvOefudc4tcM5Nl5f4/lfSn/wzZ6LbKamxlAoZT/+zbS95f5x5Xd4K5UJJb8u7dlryNrWIbqckxzMZChKnQsc0LInaV/KWb4NeyNnyln3z+utLYY+P2C7v4v+qcZ5PVj/58VVUu0F9ScFLqcVJ0uJpZqMljZS3tfStcfqKtBvUl1S841kkc9M5t0He1ry35+gr0m5QX1LxjqVUiHj6xzwl77q+/5pZYzNrLG87dEk63S+L3MOGeJ74/MyKtJ2jr0i7QX1JxHOH//W/Ac997X+NvA8Sz4K9F6VJukne9t8L4vQVaTeoL6l4x7Owsbxb3gfYV6MLnXMH5SUa9eVtzBTPDMGAAAAIHUlEQVTpK9JuUF9S8Y6llIS56Zx7Vd61+63lbSRUxzl3q192VNJnUX1F2g3qSyr+8UyGRHFyUXUKHdOwJGqr5Y2lXXShmZWX1ErSmiI+PuI8eZN2TxH3kx+r/a8dAp5rL28irE1if6mQlHj6SdooSbMkDc550azvf+UtPceLp/LbX0gV5dysIO/Uk4iSHkupcPGsIKmmvFWef0c9VvjP9/X/P9j/P/E8wdfnH3+OYjcOIp6JX19kY4JzAp6LlH0b1ZfEe1F+XSeptqQX41xXVNLnZ2FjGfnwWjrgubQcX5mb+ZwrzrljzrkPnXMrnXPfmtmZ8hK3f/hJsFTy52ayJJp3nzjnMvJZN/EcPRnbX+Zjq8vmyvs+EX2jyhpJalqI409X1Fa5UeXX+nXjbnVa2HEGtJfoPmqr5S1J14kqq+OXLU319y3V33e//GG/7gtKfOPHV+XdH6RlVFnk/iCfqhjfHyQJP0OB28NKutyP2bIc5SU2loWNp7xz/m8IeNzmH7vY///5xDPf87N6nHZ/p+D7/BHPvONZ1X8f2aGoWx1IOkve9UGf5KjPe1Ee8cxxzN/9Y+LeyLYkz88kzM2n4/xMnyFvdWKPYm93wNzM59yMqldK3kYtx+Vdo3pKzM04sUj0ObyevJ2by0SV1ZR3uuL7Cr6PWs570xVqjqY8SFGDnqIf7rw+WNKT8i50XKGoD+DyTk90hTi+u7yLpX8v7wL0O+StxByTt9vd+UU8zn7yTtEbKe+vwN9F/b9fjro/kffXjS2ShvqPLfLeSFvmNc7i8ihMPP3vnfN/gdwkb6Ui+nFljvqN5f2S/0bS/fJO51svbxW1S6pjkeJYLpC3Veyjkm7xfzZekHf9z15F3fzyVIhlYeMZp70Gin/Da+KZ9/x8Wt72xo/KO7V5uLxdH50/bysQzwK/F/3KP36jvJux3y/vd2mmpKty1OW9KEE8/efq+HPs/QR9lej5Wcif9fryThs9Lmm2//P+gLzdup2k25mbBYpnZUkfy9tsabCke+StijnluDHzqTA3/ddYkM/hK/xYNchRfo9f/rb/u3SMP+c2Kcd9Pgs7R1MesKgXUtp/4Z/4L+hLedd55HzBgb8kC3D8j+T9JSESpEjwpirqruFFOM7INz3osSKgfgdJy/yx7pe3ZXqbVH+/wvB9l3fz23ixjBfPH0l6TV7ycVDSO5I6pzoOIYjlL+T9JfgLeTsRHZK3++MUSfXi9FdiY1nYeMZpr4HiJGrEM+H8vN7/3felPz8PyLuH0AOSyhPPE5ufknrIS3QPyHt/eVPST+PU5b0ocTwf8H/Gh+SjvxI7PwsbS3krQ7PkrfhmyVt5+B9JPeL0x9yME09JZSW9LC/RPSwvCXtDeSRdJXlu+q9vhfL5uVFxEjX/uQHybih+WN6p4jMk1Ur2HDW/AQAAAABASIRlMxEAAAAAgI9EDQAAAABChkQNAAAAAEKGRA0AAAAAQoZEDQAAAABChkQNAAAAAEKGRA0AAAAAQoZEDQAAAABChkQNAAAAAEKGRA0AUOTMrKKZDTWzlWa2x8yyzOwbM1tkZgPMLC3gmKpmdsjMnJn1i9NuJ//56EeGma0zs7uD2gUAoDjgDQwAUKTMrLGk1yWdL2mppMck7ZJUS1JnSc9LukDSfTkO7SOpnKRtkm6WNDuPbl6WtEiSSTpT0k2SnpL0I0m/StJLAQDgpDHnXKrHAAAoocysgqT1khpJ+qVzbn5AnYskXeScm5ajfL2kPZJekzRJUmPn3NYcdTpJelvSvc65iVHllSRtlnS2pNrOuZ3JfF0AABQ1Tn0EABSlwZKaSHoyKEmTJOfc6oAkrY2kVpJmSZoj6ai8VbV8cc4dkPSevBW2Ric2dAAAUodEDQBQlG7wvz5bwOMGScqQ9Bfn3C5Jf5fU38wK8r4VSdD2FLBvAABSjkQNAFCUfixpX85TFvNiZuUl3SgvSTvgF8+SdI6kLnEOq2hmNcysppk1N7OpklpL+sA592khxg8AQEqQqAEAilIVSfsLeEwPSWfIS84iFknaqfinP47xn/9W0gZJt0uaL+n6AvYNAEAokKgBAIrSPkmnFfCYQfKSrh1m1tjfNbK+pDcldTOzGgHHPCvpSknXSPqtvNMdz5F0OKgDM3vPzPoWcFwAAJw0bM8PAChKGyVdZmYN83P6o5mdK+lyeZuAxDtlsa+8XSCj/ds5t9T/92Ize0fSO5KekdQrRx+lJDWX9GG+XwUAACcZK2oAgKL0F//r4HzWHygvSRsiqWfA4xPlY/dH59y78u679ksz+0mk3MzOlLdJSQVJ7/k3x74sn2MDAOCk4T5qAIAiY2YVJa2T1FBST+fcawF1LpR0sbzVr+2S9jrnWsRpb5Sk0ZLaOedWx7uPml+3sbx7qa1wznWOKu8labhzrm2hXyAAAEWEFTUAQJFxzh2U1FXSNkl/NbM3zGy4mQ00s/vMbLGk1ZLqSbpKUl39sAoXJPLcoHz0/ZmkVyRdYWaXRj3VRl7yCABAaJGoAQCKlJ8wtZY0TFIlSQ/K2/zjHknHJfX3yyLJV+CNsf22Nsq7dq2XmVXIR/eP+H2MjSojUQMAhB6nPgIATilmtkvSNc65D1I9FgAA4mFFDQBwyvB3fDwj1eMAACAREjUAwCnDOXdc0kRJb/o7PjZO9ZgAAAjCqY8AAAAAEDKsqAEAAABAyJCoAQAAAEDIkKgBAAAAQMiQqAEAAABAyJCoAQAAAEDIkKgBAAAAQMiQqAEAAABAyJCoAQAAAEDI/B/RGjaRrmKGzgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "font = {'family' : 'sans-serif',\n",
    "        'weight' : 'normal',\n",
    "        'size'   : 18}\n",
    "\n",
    "matplotlib.rc('font', **font)\n",
    "\n",
    "def simpleaxis(ax):\n",
    "    ax.spines['top'].set_visible(False)\n",
    "    ax.spines['bottom'].set_visible(False)\n",
    "    ax.spines['right'].set_visible(False)\n",
    "    ax.spines['left'].set_visible(False)\n",
    "\n",
    "c_ = sns.color_palette('Greens')\n",
    "\n",
    "f,ax = plt.subplots(1,1,figsize=(15,5))\n",
    "\n",
    "## CAR(week 0)\n",
    "sns.distplot(dfreport['median CAR0'],color = c_[5],ax=ax, kde=False, norm_hist=True, label='region-wide week 0')\n",
    "mu = np.median(dfreport['median CAR0'].values)\n",
    "ax.axvline(mu, color='g',label='median week 0')\n",
    "print(mu)\n",
    "print(dfreport['median CAR0'].quantile(q=[0.025,0.5,0.975]))\n",
    "\n",
    "\n",
    "##CARc\n",
    "sns.distplot(dfreport['median CARc'],color = c_[2],ax=ax, kde=False, norm_hist=True,label='region-wide April 15th, 2020')\n",
    "mu = np.median(dfreport['median CARc'].values)\n",
    "ax.axvline(mu, color='g',alpha=0.5,label='median April 15th, 2020')\n",
    "print(mu)\n",
    "print(dfreport['median CARc'].quantile(q=[0.025,0.5,0.975]))\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "ax.set_xlim((0.001,0.6))\n",
    "plt.legend()\n",
    "\n",
    "# ax.set_title(label,fontsize=fs)\n",
    "# ax.axvline(mu)\n",
    "# ax.text(mu,20,str(np.round(mu,3)))\n",
    "ax.get_yaxis().set_visible(False)\n",
    "ax.set_xlabel('')\n",
    "plt.subplots_adjust(hspace=0.5,top=0.8,bottom=0.2)\n",
    "\n",
    "def calcfold(x):\n",
    "    return 1/x\n",
    "def calcfold2(x):\n",
    "    return x\n",
    "# ax.secondary_xaxis('top', functions=(calcfold,calcfold2))\n",
    "x = np.array([0.05]+list(np.arange(0.1,1.1,0.1)))\n",
    "ax.set_xticks(x)\n",
    "ax.set_xlabel(r'CAR$_{t}$')\n",
    "ax2 = ax.twiny()\n",
    "ax2.set_xticks(x)\n",
    "ax2.set_xticklabels(np.round(1/x,1))\n",
    "ax2.set_xlabel('fold difference')\n",
    "simpleaxis(plt.gca())\n",
    "plt.savefig(productpath + 'carpastpresent.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
